

%==============================================================
% Driver program to calculate the V_coinc and V_coherent for a 
% network of gravitational wave detectors.
%==============================================================

clear all;

                   
  % make a source sky grid with good resolution:
  angularResolution = 0.05;
  skyPositions = sinusoidalMap(angularResolution);

  % extract the spherical angles associated with the sky grid:
  thetas = skyPositions(:,1);
  phis = skyPositions(:,2);

  detNames = {'LHO','LLO','Virgo','INDIGO'};
  nDetectors = length(detNames);

  psi = 0.0;
  iota = 0.0;
  mTot = 2.8; %double NS binary
  indigoDetPsi = 0.0; %degrees

  [rmax2, rmax3] = Rmax_coincident_comb(iota, psi, phis, thetas, detNames, indigoDetPsi, mTot);

  rmax2_coinc_HLV  = max( [rmax2.HL; rmax2.HV; rmax2.LV] );
  rmax2_coinc_HLVI = max( [rmax2.HL; rmax2.HV; rmax2.HI; rmax2.LV; rmax2.LI; rmax2.VI] );

  coinc2_f_ratio = (rmax2_coinc_HLVI - rmax2_coinc_HLV)./rmax2_coinc_HLV;

  rmax3_coinc_HLV  = rmax3.HLV;
  rmax3_coinc_HLVI = max( [rmax3.HLV; rmax3.HLI; rmax3.HVI; rmax3.LVI]);
  
  coinc3_f_ratio = (rmax3_coinc_HLVI - rmax3_coinc_HLV)./ rmax3_coinc_HLV;

  % Horizon volume in Mpc^3
  dOmega = skyPositions(:,4);

  coinc2_vol.HLV = (1/3) * sum ((rmax2_coinc_HLV).^3 * dOmega);
  coinc2_vol.HLVI = (1/3) * sum ((rmax2_coinc_HLVI).^3 * dOmega)

  coinc3_vol.HLV = (1/3) * sum ((rmax3_coinc_HLV).^3 * dOmega);
  coinc3_vol.HLVI = (1/3) * sum ((rmax3_coinc_HLVI).^3 * dOmega)

  % Effective spherical radius
  coinc2_Reff.HLV = (3/(4*pi)) * (coinc2_vol.HLV)^(1/3);
  coinc2_Reff.HLVI = (3/(4*pi)) * (coinc2_vol.HLVI)^(1/3)

  coinc3_Reff.HLV = (3/(4*pi)) * (coinc3_vol.HLV)^(1/3);
  coinc3_Reff.HLVI = (3/(4*pi)) * (coinc3_vol.HLVI)^(1/3)

  %========== Plotting code ============
  try  
      Plot_Rmax_Coincident;
  catch
      fprintf(2, 'Could not execute Plot_Rmax_Coincident\n');
  end

  %================================================================================================
  %================================================================================================
  
  k = 1;
  indigoDetPsi = 0; %degrees
  % Vary the detector's arm geometry and plot the change in Rmax
  while indigoDetPsi <= 90
       
       clear rmax2 rmax3 rmax2_coinc_HLVI rmax3_coinc_HLVI  ;  
       [rmax2, rmax3] = Rmax_coincident_comb(iota, psi, phis, thetas, detNames, indigoDetPsi, mTot);
        
       rmax2_coinc_HLVI = max( [rmax2.HL; rmax2.HV; rmax2.HI; rmax2.LV; rmax2.LI; rmax2.VI] );
       rmax3_coinc_HLVI = max( [rmax3.HLV; rmax3.HLI; rmax3.HVI; rmax3.LVI]);

       mean_rmax2_HLVI(k) = mean( rmax2_coinc_HLVI );
       mean_rmax3_HLVI(k) = mean( rmax3_coinc_HLVI );

       k = k+1;
       indigoDetPsi = indigoDetPsi + 5;
  end
  figure(167); clf;
  plot([0:5:90], mean_rmax2_HLVI, 'b-', 'LineWidth', 2); 
  hold on;
  plot([0:5:90], mean_rmax3_HLVI, 'm-', 'LineWidth', 2); 
  hold off;
  title('Average Coincident Distance vs. IndIGO Arm Geometry -- aaaa' , 'FontSize', 12);
  xlabel('IndIGO Arm Geometry / Deg', 'FontSize', 14);
  ylabel('Avg. Coinc. Dist. / Mpc', 'FontSize', 14);
  axis tight; box on;
  set(gca, 'FontSize', 12);
  legend('2-IFO', '3-IFO');
  print('-dpng', 'Coincident_VaryDetPsi_aaaa.png');

  %================================================================================================
  %================================================================================================
 
  k = 1;
  indigoDetPsi = 0; %degrees
  mTot = 2.0; %solar masses 
  clear mean_rmax2_HLVI mean_rmax3_HLVI;
  % Vary the masses and plot the change in Rmax
  while mTot <= 200
       
       clear rmax2 rmax3 rmax2_coinc_HLV rmax2_coinc_HLVI rmax3_coinc_HLV rmax3_coinc_HLVI  ;  
       [rmax2, rmax3] = Rmax_coincident_comb(iota, psi, phis, thetas, detNames, indigoDetPsi, mTot);
        
       rmax2_coinc_HLV  = max( [rmax2.HL; rmax2.HV; rmax2.LV] );
       rmax2_coinc_HLVI = max( [rmax2.HL; rmax2.HV; rmax2.HI; rmax2.LV; rmax2.LI; rmax2.VI] );
       rmax3_coinc_HLV  = rmax3.HLV;
       rmax3_coinc_HLVI = max( [rmax3.HLV; rmax3.HLI; rmax3.HVI; rmax3.LVI]);

       mean_rmax2_HLV(k) = mean( rmax2_coinc_HLV );
       mean_rmax2_HLVI(k) = mean( rmax2_coinc_HLVI );
       mean_rmax3_HLV(k) = mean( rmax3_coinc_HLV );
       mean_rmax3_HLVI(k) = mean( rmax3_coinc_HLVI );

       k = k+1;
       mTot = mTot + 10;
  end
  figure(267); clf;
  plot([2:10:200], mean_rmax2_HLVI, 'b-', 'LineWidth', 2); 
  hold on;
  plot([2:10:200], mean_rmax2_HLV, 'b-o', 'LineWidth', 2); 
  plot([2:10:200], mean_rmax3_HLVI, 'm-', 'LineWidth', 2); 
  plot([2:10:200], mean_rmax3_HLV, 'm-o', 'LineWidth', 2); 
  hold off;
  title('Average Coincident Distance vs. Total Mass of Binary System -- aaaa', 'FontSize', 12);
  xlabel('Mass_{Total}/ M_{sun}', 'FontSize', 14);
  ylabel('Average Coincident Distance / Mpc', 'FontSize', 14);
  axis tight; box on;
  set(gca, 'FontSize', 12);
  legend('2-IFO (HLVI)', '2-IFO (HLV)', '3-IFO (HLVI)', '3-IFO (HLV)');
  print('-dpng', 'Coincident_VaryMtot_aaaa.png');
 
   
    
    
   



                    
                    
     
        







